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Abstract. We offer an alternative viewpoint on Dyson’s original paper 
regarding the application of Brownian motion to random matrix theory (RMT). 
In particular we show how one may use the same approach in order to study 
the stochastic motion in the space of matrix traces t n = rather than 

the eigenvalues X u . In complete analogy with Dyson we obtain a Fokker-Planck 
equation that exhibits a stationary solution corresponding to the joint probability 
density function in the space t = (t \,..., £ n ), which can in turn be related to the 
eigenvalues A = (Ai,..., Ajv). As a consequence two interesting combinatorial 
identities emerge, which are proved algebraically in the appendix. We also offer 
a number of comments on this version of Dyson’s theory and discuss its potential 
advantages. 


1. Introduction 

In his seminal 1962 paper, A Brownian-Motion Model for the Eigenvalues of a Random 
Matrix [I], F. Dyson provided a conceptually novel and practical approach to the 
theory of random matrices, paving the way for many interesting developments (see 
e.g. mg and references cited therein.) In it he explains how to introduce a dynamical 
approach to the theory of random matrices and the traditional Gaussian ensembles in 
particular. We briefly recapitulate the results here in this introductory section. 

Consider a self adjoint matrix M of size N x N, whose entries are of the form 
Mij = Mjj :a e a . The coefficients Mij ;a being real parameters and e a are the 

units of the three potential algebras: real (/3 = 1), complex (/3 = 2) and real- 
quaternion (/3 = 4), satisfying eg = 1 and ej = -1 V « > 0. Choosing the real 
coefficients Mi ra independently from a Gaussian distribution with zero mean and 
variance E = (1 + fy,)/(2/?) we obtain the Gaussian orthogonal, unitary and 
symplectic ensembles (GOE, GUE and GSE) for /? = 1,2 and 4 respectively. Thus the 
probability distribution for the matrix M may be neatly summarised in the following 
form 

P(Af) = nf ) e -§ trMMt 

with a normalization constant. 


( 1 . 1 ) 
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Crucially, Dyson realised that the above distribution can be identified as the 
stationary distribution of a Brownian particle in N + /3N(N — l)/2 dimensions. More 
precisely, this means that each independent element Mij. a , 1 < i < j < N undergoes a 
ID Ornstein-Uhlenbeck process, so that in the (fictitious) time s the motion of Mij- a 
is completely determined by the following moments: 

E {SMij-c) = - 5s (1.2) 

E(JM£.J = i(l+ *«)&»■ (1.3) 

The latter implies that E(|(5Myj 2 ) = (1 + (2//3— 1 )6ij)5s (since the diagonal elements 
Mu are always real). 

Importantly, this stochastic motion is invariant under unitary transformations, 
meaning the eigenvectors do not play any role in the corresponding motion induced in 
the N dimensional space of eigenvalues A = (Ai, • ■ •, A tv). Therefore one may choose 
a representation in which M is diagonal, leading to a perturbation of the eigenvalue 
A^ due to a small change in the matrix SM of 

^ = + (1.4) 

, An - Ap 

Obtaining the first two moments of the evolution then follows directly from the 
expressions (E2D and (E3D, given by 


E(JA M ) = P M (A)<5s 
E(5A£) = Its. 


E 


1 


A,-A 


v 




6s 


(1.5) 

( 1 . 6 ) 


Using these two moments, one obtains a Fokker-Planck equation that describes how 
the joint probability distribution function (JPDF) P(A; s) evolves in time, given some 
specific initial distribution P(A;0); 


dP 


N 

E 


-^-(P,(A)P(A; S )) + /3- l52P(A;S) 


dX 2 , 


(1.7) 


The real advantage, and one might add elegance, of this approach is expressed in the 
above equation. In general it is not known how to obtain P(A; s) for arbitrary initial 
conditions and times s. However, since we are interested in the stationary distribution, 
we can reduce the complexity by setting the LHS equal to zero, at which point one 
solves the equation easily to obtain 


p(A) = Cp n IV - A, 


exp 


li<v 



( 1 . 8 ) 


with Cp N ' > a normalisation constant (see e.g. Chapter 3 of [1]). Moreover, since we 
know that the underlying motion (11.21) and m in the space of matrices leads to 
the probability distribution m, the expression (II.SI) must be the unique stationary 
distribution for the process (ffUD and (11.61) and is therefore the JPDF of the eigenvalues 
in the appropriate Gaussian ensembles. 

The key component of (11.81) is the Vandermonde determinant ru, i-v - -M> 
which is responsible for the apparent repulsion of neighbouring eigenvalues. This 
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factor emerges as the Jacobian of the transformation from (11.11) to (11.81) . However, as 
Dyson highlights, the above approach offers a new insight into its appearance - as it is 
nothing more than the effect coming from the second order term in the perturbation 
formula (El. 

Recently, the authors have adapted the above approach to investigate the spectral 
statistics of Bernoulli matrices |17j (matrices in which the elements come from the 
set {±1}). In this instance higher terms in the perturbation formula had to be 
accounted for, which meant assumptions regarding the delocalisation of eigenvectors 
were required. This inevitably led to the following question - can Dyson’s Brownian 
motion model be used without the requirement of the perturbation formula (11.41) ? 

In this article we demonstrate that the answer is indeed positive. To achieve 
this we start from a slightly different viewpoint to Dyson: Rather than following the 
evolution P( A; s) of the eigenvalues directly, we instead follow Q(t; s ) - the JPDF of 
the iV-dimensional vector of traces t = (ti, • • •, tjv), where t k = = tr M k . 

Performing a transformation of variables then allows us to recover the stationary 
solution (11.81) expressed in terms of the t variables. 

To the best of our knowledge, the distribution of the traces (or spectral moments) 
has not been extensively studied, although there are exceptions for both the Gaussian 
and circular ensembles (see e.g. EHI2 and references therein). We therefore find it 
worthwhile to pursue this direction, not only as it sheds new light on Dyson’s approach, 
but because it may offer different perspectives on such trace distributions. In addition, 
our method has led to the discovery of two identities (see Proposition |T]) that relate 
the traces t n with n > N to those with n < N. We are unaware of the existence of 
similar identities in the literature and a direct proof of their validity has kindly been 
supplied by D. Zagier in | Appendix A[ 

It has also been brought to our attentioiQ that a similar philosophy has also been 
undertaken by Bakry and Zani m- Rather than looking at the traces t they follow 
the motion of the secular coefficients (given by c k in Section O . Their motivation 
comes from wanting to generalise the probability density functions to Gaussian random 
matrices with Clifford algebras (rather than real, complex or quaternion entries) and 
they too note that such approaches have not been utilised before. 

The paper is organised as follows: In Section[2]we introduce the basic concepts to 
be discussed, provide some useful relations and outline the identities mentioned above. 
We also provide explicit formulae for the stationary distribution Qp(t) for j3 = 1,2,4 
and arbitrary dimension N. In Section [3] we derive the Fokker-Planck equation for 
Q(t\ s) and give an example of its form in two dimensions in Section r3.ll In Section [3.21 
we analyse the equation in ?V-dimensions and show how the aforementioned identities 
arise from considering the stationary solution Qp(t). Section EP1 is used briefly to 
explain how the mean spectral moments also arise naturally in this context. Finally 
in Section [4] and Section [5] we provide an application of this method to Bernoulli 
ensembles and discuss the potential advantages of the whole approach. 

2. Definitions and useful relations 

The first essential feature to outline is the relationship between the spectral and 
trace distribution functions P(A;s) and Q(t;s). The elements of the Jacobian of 

| For which we would like to thank P. Forrester. 
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the transformation are given by = nWf , which means that 

dt 


P(X;s) = 


d\ 


Q{t; s) = N\ det(V) Q(t; s ). 


( 2 . 1 ) 


Here V is the familiar Vandermonde matrix 
/l 1 ••• 


V = 


Ai 

A? 


^2 

X 2 

A 2 


1 \ 

A N 


\2 

a n 


V Ai' 


N-l 


A 


N-l 


A 


N-l 

N 


( 2 . 2 ) 


and so det(V) = n,,<„ | A M — |, as seen in m ■ The mapping A —> t is one- 
to-one as long as the Jacobian does not vanish, hence we must restrict the spectral 
variables to an ordered sector, e.g. Ai < A 2 < • • • < Xn- In order to obtain an 
expression for Q(t; s ) we need to write P(A; s), and thus the Vandermonde determinant 
det(V), in terms of the traces t. Fortunately this is relatively straightforward, since 
G(t) = det(V) = i/det(W T ), with 



( to 

tl 

t/v-l \ 


tl 

t2 

tN 

VV T = 

t2 

£3 

tN+1 


\ tN-1 

tN 

■ ■ ■ t 2 N-2 / 


and to = N (see PUCES] for example for uses of this identity in other contexts). 
At this point G{t) is expressed entirely in terms of the traces, as desired, however 
this includes traces of higher degree than N , which are themselves functions of the 
traces t n , 1 < n < N. The expressions for tN+r in terms of the first N t n , whilst 
complicated, can be written down explicitly. They originate from the characteristic 
polynomial $(V) := det(XJ — M) = Y^k=o c kX N ~ k , with Co = 1. For any eigenvalue 
X v we have $(A l/ ) = 0 and thus it follows 


tN+r — 


N 


^ CktN+r-k 


Lfc=l 


(2.3) 


Newton’s identities give the coefficients Ck in terms of the t n ,n < k via the determinant 


Cfc = 


k\ 


ti 

ti 

tk-X 

tk 


1 

tl 

tk -2 
tk -1 


tl 

t2 


0 

0 

k-1 

tl 


(2.4) 


Therefore, using a combination of relations (USD and not one may write G(t ) 
explicitly in terms of the first N traces t. Clearly A = G(t) 2 is nothing but the 
discriminant of $(Af) expressed as a function of t. 

Using the transformation EH) and the stationary distribution for the eigenvalues 
(11.81) we can obtain the JPDF for the traces in the three canonical ensembles 

fs.\ n^ N '>n(+\P- 1 @ ^ 
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Xjv(t) is an indicator function for the domain T C which is the support for Qp(t). 
In contrast to the spectrum, which is defined over the entire space R^, the trace 
parameters are restricted to the domain T. This is because the traces are sums of 
powers of real variables, which must satisfy certain consistency relations: The inverse 
mapping t —> A should yield real spectra. For example, in 2 dimensions we have 
2t 2 — t\ = (A 2 — A 1) 2 > 0. Hence, T = {(ti, ^ 2 ) G R 2 : 2t 2 — t\ > 0}. In higher 
dimensions it becomes increasingly more difficult to write an explicit definition of 7~, 
other than stating that it is the image of the mapping A —> t. It should be emphasized, 
however, that T is independent of the ensemble under consideration - one may consider 
matrices with non-Gaussian elements, or even correlated elements, and T will remain 
the same. 

We would also like to highlight that the GOE distribution takes a very simple 
form in this space, i.e. Q\(t) = c[ N ^ exp (—^ 2 )- At first sight it might seem strange 
that the JPDF for all the traces depends only on one parameter ; however, as alluded 
to above, one must pay very close attention to the domain of integration T. This is 
exemplified in Section 13.11 in which we calculate expectations values and marginal 
probabilities. 

In the following section we shall derive the Fokker-Planck equation for Q(t;s). 
Its stationary solution is known and given explicitly in (12.51) . As will be shown below, 
by substituting this solution into the stationary Fokker-Planck equations we obtain 
two identities which are summarized in the following proposition. 

Proposition 1. For n > 0 we have 
N 

2 171 + (rc + l )tn (2.6) 

™ — 1 m i,j> 0 

i+j = n 

= ^2 Utj - (n + 1) . (2.7) 

i, j >0 
i+j=n 

As mentioned in the introduction, we are unaware of such identities arising before 
in RMT or any other context and a direct algebraic proof is given by D. Zagier in 
|Appendix A| 

3. The Fokker-Planck equation 

The main reason for studying Dyson’s Brownian motion in the space of traces is that 
the Fokker-Planck equation for Q(t; s ) can be derived directly, avoiding the use of 
perturbation theory USD. The expectation values of the components of t due to 
an incremental changes in the matrices will be evaluated directly from the matrix 
elements statistics. Once Q{t\ s ) has been computed, one can then transform back to 
the spectral representation in order to deduce the eigenvalue statistics though P( A, s). 

We begin by expressing the change in the n th trace via the change in the matrix 
6M, up to second order (since higher terms in SM will be of orders Ss 2 or greater 
after taking the expectation) 

5t n = [tr((M + SM) n ) - tr(AT 1 )] 

n —2 

= ntr(M n ~ 1 6M) + -J2 tr(M x dMM n - x ~ 2 SM) + ... . (3.1) 

2 *=o 


N 


2 ^ 3G_ 

— 2^ m tn+m 


m— 1 


dt n 





Dyson’s Brownian-motion model for random matrix theory - revisited 


6 


The simplest way to compute E(dt„) and E(4 n c)f m ) is to invoke the invariance of the 
stochastic motion under unitary transformations. We are then free to write the initial 
matrix Min a diagonal representation of eigenvalues, i.e. M, t] = XiSij. Using this and 
the expressions and m we find 


E(5t„) 


-nt n 8s+^Y,Y. 


2 4zE (5M jk 5M H ) 


x=0 ijkl 
n—2 


ll \ 

?T'tn H"~ 77 / tx^n— 2—a 


x—0 


— -{n-l)t n _ 2 


8s, 


(3.2) 


where we have used that ’E(\8Mij\ 2 ) = (1 + (2//? — l)Sij)Ss. In particular this means 
for n = 1 and 2 that we have E(4i) = — t±8s and E(4 2 ) = (—2t 2 + to + (2//3 — l)io)<5s. 

For the second order moments, since again we need terms proportional to 8s and 
no more, we only require the first term in (13.11) . Therefore, for n,m = 1,... ,7V, we 
get 


E (Stjtm) = nm'52%-%Xk- ls ME(6M ji 5M lk ) 

ijkl 


= nmY J K~ 1 K~ 1 ^ M ^M kk ) 

ik 


2 nm 

~r 


t 


n+m 


-2 8s, 


(3.3) 


where we have used E(5M i 4M fc j.) = j8ik8s. Note that in the above equations, and in 
the following, one should remember that the independent parameters in the present 
theory are the components t which consist of the first N traces. Whenever there 
appears t x with x > N, it should be considered as a function of the independent 
parameters as explained in the previous section. Similarly, one must substitute to = N. 

We are now in a position to obtain our Fokker-Planck equation for determining the 
probability distribution Qp(t;s) of the traces. For simplicity we write (13.21) and (13.31) 
in the form R^ = E(8t n )/8s and = E (8t n 8t m )/8s, so that (see for instance [23 ] ) 

dQp ^ djR^Qp) l^ d 2 {R { n^Qp) ( 

ds ^ dtn + 2 ^ dt n dt m ' 1 ' ’ 

n n,m 


Just as Pp(\), given in (11.81) . is the stationary solution to the Fokker-Planck equation 
CD for the eigenvalues, so we would like to verify Qp(t), given in (12.51) . is the 
stationary solution of (13.41) above. For this to be the case, Qp(t) must therefore 
satisfy the following N simultaneous equations 

= ir 1 . V !<»<». (3.5) 


These will be discussed shortly for arbitrary matrix dimension N but prior to this we 
outline, for illustrative purposes, the scenario for N = 2. 


3.1. Example: 2x2 Gaussian ensembles 

The N = 2 case offers the particular advantage that the expressions (13.21) and (13.31) do 
not contain traces larger than tjv (i.e. f 2 in this case), which is not true for N > 2. In 
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order to satisfy (13.41) Q = Qp(ti,t 2 ) must be a solution of the simultaneous equations 
3751) . which in 2 dimensions are given by 


0 — tiQ + 


1 

1 L 


(t 0 dQ) d(iiQ) 


0 = ( 2fo — ti — 


dt\ 

(2-/3) 


+ 2 - 


dt 2 


to ) Q + 


d(tiQ) , d(t 2 Q) 


dti 


+ 2 - 


dtn 


One may verify by substitution that the solution is, including the normalisation 
constant presented in (11.81) . 


Qp{t\,t 2 ) 


2 C t ] ( 2 * 2 -*?)' 


0*2 


(3.6) 


Written in terms of the eigenvalues, using G(i ) 2 = (2f 2 — t 2 ) = (A2 — A1) 2 , this yields 

P/3(Ai, A 2 ) = G^ 2 ) |A 2 - A 1 |Q(A 1 , A 2 ) = cf ] \\ 2 - A 1 |^- £ ^ tMl , 

which is the expected result for the JPDF. 

From (13.61) we can immediately calculate the marginal probability distributions 
for the traces. Importantly, the limits of integration are defined by the domain T. For 
2 dimensions this was outlined in Section [2] 

qp(ti) = / dt 2 Qp(t 1: t 2 ) = oCl 2) r73e _/3t i /4 (3.7) 

Jtl/2 2 

/ v/2*I -1 

dfi Qp{ti,t 2 ) = -Cf ) spt^ /2 e~ pt2/2 , (3.8) 

-V2t7 2 

where (G^) -1 = 4-^/tF, 7 r, 37r/8, rp = 2, and sp = 2 3 / 2 ,7r, 37r/2 for 

/3 = 1, 2,4 respectively. The expected value of t 2 is therefore (t 2 ) = / 0 °° dt 2 t 2 qp(t 2 ) = 
3, 2, 3/2 in the three cases. 


Stationary solution 


Finding the stationary solution in N dimensions requires solving the N simultaneous 
equations given by (13.51) . Therefore, substituting in the expressions (13.21) and (13.31) we 
get for each n 


Tl \ t 

Tltn + 77 / v t x tn—2—a 


x—0 


N 


2 — 3 n . \ ii v—> 

—g- 2 ~ l)tn ~ 2 1 ^ = 7 S 


m a(t n+m _ 2 Q^) ^ (3g) 


The derivative in the RHS can be expanded using the chain rule to obtain 
9(tn+ m - 2 G /3_1 e _/3t2 / 2 ) 


dtr 


dt n+ m-2 , + (/3-l)<9G /3^ x 

dt m +tn+m ~ 2 G dt n 2 n+m ~ 22m 


Qp- 


Therefore, after some algebra in which we divide through by a factor nQp/(2(3) and 
cancel the term involving — nt n on both sides, we arrive at the following relationship 
between the traces 

n —2 N r 

P ^2 t x t n - 2 - x +{ 2 -P)(n-l)t n -2 = 2 ^2 rr Ln+m ~ 2 


x—0 


dtr, 


+ tn+m—2 


(/3 - 1) dG 


G dtr, 


.(3.10) 
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In the particular case /? = 1 there is no dependence on the Vandermonde determinant 
G(t) and we get 


N 




dt 


n—2 


n+ra—2 
dim 


— ^ ( txtn—2—x T + l)^n—2- 


(3-11) 

m= 1 ~ " L x=0 

Replacing n — 2 by n thus gives the identity (12.61) . If we then rearrange (13.101) in terms 
of /3 we find 

/ n—2 N 


/3 E txtn ~ 2 ~x ~ ( n - 1)^-2 - 2 y ^ mt 


n+m—2 


^a:=0 


1 dG 
GdtZ 


= 2 



dt n+ m -2 + 1 dG 

dt m tn+m ~ 2 Gdt m 


(n 



(3.12) 


The above must be fulfilled simultaneously for both /3 = 2,4, which only occurs if the 
expressions in large brackets on the two sides of (13.121) vanish. Therefore, by using 
the substitution (13.111) we arrive at the second identity (12.71) 


o A 1 dG 

2 / Ultn+m—2 


m =1 


Gdt n 


n—2 

E 

x=0 


— / J txtn—2—x {p> l)^n—2 


where again we must replace n — 2 by n. Since we know that the expression m 
must be our stationary solution the method above constitutes a proof of the identities 
(12.61) and (12.71) . However, a direct proof of these is given by D. Zagier in |Appcndix A[ 
which therefore implies that (12.51) must be our stationary JPDF, without the need for 
any transformation of variables. 


3.3. The mean values ( t n ) 


Computations of expected values of any function of t involve integrating over the 
domain XTv(^)> which is not explicitly defined for any N > 2. However, one can use a 
simple heuristic reasoning in order to identify the mean values ( t n ) as the coordinates 
of the vector t for which the drift force (13.21) vanishes, i.e. 


(tn) ~ 2 ^ 

X—0 


1 O _ ft 

^(t x )(t„- 2 - x ) + —+(n - l)<i„_ a ). 


2/3 


(3.13) 


It is natural, and customary, to scale the matrices M by 1 /y/N and the resulting 
traces by 1/N, so that we may define r n = N~^^ 1 t n - Thus 

i n ~ 2 o _ a 

(Tn) = 2 ^2(Tx)(t„.-2-x) + -^-( n - 1) (Tn—2) (3.14) 

Z x=0 P 

If we take /3 = 2, with initial conditions To = 1 and T\ = 0, then (13.141) implies 
that (T 2 fc+i) = 0 and +&} = where Ck are the Catalan numbers. This is the 
well known result obtained by computing the moments using the semi-circle spectral 
distribution function (see e.g. sil[9][l8]). For other /3 the last term is of order 0+) 
smaller than the rest and thus its effect vanishes in the limit of large N. This is 
consistent with the fact that the spectral distribution of the three canonical ensembles 
converge to the semi-circle distribution for TV — > oo. Moreover for N = 2, (13.131) 
returns +} = 3,2,3/2 for /? = 1, 2,4 respectively, which is exactly the result obtained 
in Section l3Tl 
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Recently, the authors have used a discrete analogue of Dyson’s Brownian motion model 
to investigate the spectral statistics of Bernoulli ensembles mi. Here we provide a 
brief illustration of how this can be adapted to the traces setting and discuss why this 
offers certain advantages. Our Bernoulli ensemble 21 at is given by the set of N x N 
symmetric matrices with 0 on the diagonal and off-diagonal entries chosen randomly 
and independently from the set {±a} with equal probability (in the following we shall 
choose, without loss of generality, a = \j\[2 in order to match the variance of the 
GOE defined in Section [[]). The spectral properties of 21 at were first analysed by E. 
Wigner in 1955, who showed the empirical spectral density converges to the semicircle 
distribution in the limit of large N [18] . Recent works have gone much further, 
establishing that local eigenvalue correlations do indeed converge to the corresponding 
Gaussian expressions as N increases [T9l - f22l . 

In [f7] the random walk is defined on 21 at such that at each single time-step, one 
of the djv = \N(N — 1) off-diagonal matrix entries B pq is chosen at random and its 
sign is flipped (together with B qp ). This leads to a change in the matrix B of 


SB pq = —2B pq [\p) (q\ + |g)(p|], 


(4.1) 


where | p) is a vector whose elements are all zero but for 1 in the position p, and (p| is 
its transposed. This perturbation in turn induces a change in the eigenvalue X p of 


SX„ = (p\SB pq \p) 


E 


\{v\SB pq \p)f 
X p — A„ 


(4.2) 


in a similar manner to (11.41) . In order to construct the coefficients in the Fokker-Planck 
equation one has to average JA^ over the entire neighbourhood of matrices that can 
be reached in a single step. In particular, F,((p\6B\p)) = —2X p /dN and 


n\(v\5B\p)\ 2 ) = 


—E 

djv ■“ 

p<q 


\{v\6B pq \p)\ 2 = 


dn 


1 + dvu — 2 


N 

E 

p=i 




(4.3) 


Here, in contrast to m, there is an additional term Y2 P =i that cannot be written 
purely in terms of the eigenvalues, meaning the motion is not autonomous. 

Collating the above expressions allows one to derive a Fokker-Planck equation 
which describes the motion of a suitable observable, up to an error that depends on 
N. This error comes from a combinations of factors such as higher moments E(<5A]() 
and higher terms in the perturbation formula (14.21) . This is because, ultimately, our 
process in discrete and, unlike Dyson’s Brownian motion, one cannot assume that the 
change of the matrix due to a single step can be made arbitrarily small. These errors, 
together with the correction to the second moment from the additional term in (14.31) . 
all depend on the eigenvectors and can only be assumed to become negligible in the 
large N limit if they are sufficiently delocalised. For the present ensemble, it has been 
proved this is correct with high probability (see [T9ll22l and references therein) but for 
Bernoulli ensembles with correlated matrix entries there are no rigorous results thus 
far in this direction. Moreover, perturbation theory only converges when |(/x|(5i?|p)| is 
small relative to (A^ — A m ±i|. In ensembles such as random regular graphs, this is not 
the case, even though the eigenvectors are delocalised, due to the growth rate (or lack 
thereof) of the mean level spacing. These observations therefore motivate the search 
for another approach. 
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In complete analogy to Section [3] we can also study the random walk in the 
space of traces. In fact we shall find it more amenable to use the rescaled traces 
T n = 7V _n ' 2 _ 1 t„, as used in Section EP1 In this basis all the variables are 0(1) in N 
and thus it becomes transparent as to which terms can be neglected. To facilitate this 
transition let us therefore scale the original matrices by B = B/y/N. Applying this 
to m we have 


Sr n = -(Tr (B + SB) n - Tr (B n )) 


1 

N 


n—2 


Tr(B n ~ 1 SB) + -^ Ti -{B x 5BB n ~ 2 ~ x 8B) + ... 


x—0 


■ (4-4) 


Although we shall eventually seek to neglect those higher terms, as in Section [3] 
the whole expansion is finite for fixed n and thus exact. It means this formalism offers a 
distinct advantage over the perturbation formula (Id. 21 ) . which has no such guarantees. 
Moreover, in this way the change in the variables can be expressed directly in terms 
of the matrix elements, which is not the case for the eigenvalue representation, since 
it relies on the appearance of the eigenvectors. 

Proceeding in a similar manner, the expected change of ST n in one time step may 
be calculated as follow^ 


E(Tr(.B” - 1 AB)) = - — £ B pq TV (B^[\p)(q\ + |«)<p|]) 


p<q 


_1 Vr B n ~ l = —— N~ n / 2 t —-—Nt 


p<q 


du 


dw 


(4.5) 


and 


E (Tr(B x 6BB n ~ 2 ~ x 6B)) = — £ B 2 pq Tr(B x [\p)(q\ + \q)(p\]B n ~^[\p)(q\ + |g><p|]) 


p<q 


d N N 


n—2—x I r>x on— 2—x 
qp ' D pp D qq 


p=£q 


2 

djv 


N 




n—2—x I r>x nn-2-x\ _ \ r tdx o' 

qp ' n pp n qq ) jy n pp U ] 


n— 

pp 


where 


= 3 — [t „-2 + NT x T n — 2 —x - 2 £(x, n - 2 - x)] 
(in 


= A- 


The most striking difference between (14.61) and the Gaussian equivalent (13.21) is the 
appearance of the term f(x,n — 2 — x), which cannot be expressed in terms of the 


or 

jp °pp 


Tf Ea B qq - B PP We See tll£lt C(A S) 


variables t. Writing T s T r —C,{r, s) = 
is very close to T s T r if the diagonal elements B ® 
whole diagonal B qq . Using Wigner’s combinatorial method of counting Dyck paths 
(see e.g. EES) one can show that by averaging over 03 jv we have for fixed r and s 


are close to their average over the 


§ We use the convention that B™ q denotes the p,q -th element of the matrix B n and ( Bp q ) n is the 
matrix element B pq raised to the n-th power. 


2-x 

(4.6) 
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that (r r r s — C(r,s))<s N tends to 0 as N —> oo. Moreover, using the same technique 
one finds Var^(T r r s — £(r, s)) = 0(N~ 2 ). Hence with high probability £(r, s) is 
0(1). This shows that (Id. 61) is dominated by the term Nt x t u _ 2 -x ■ In addition 
we can also estimate those higher terms in the expectation E(<5r„) coming from the 
expansion (14.dl) . For example, we have E(Tr(<5B 3 .B ra_3 )) = lV _ra / 2 E(Tr(5i? 3 H ra_3 )) = 
4N~ n / 2 'E(Tr(SBB n ~ 3 )) = —87V _n / 2 t„_ 2 /djv = — 8T n _ 2 /dN- This again is of an order 
in N less than the dominant term in (14.61) . Therefore, in the large N limit we find 
that E(<5r n )/<5s (taking Ss = 2/cbv) tends to the expression (13.21) calculated for the 


GOE. 


Similarly, for the second moment we find 


E(ST n ST m ) = ——Y { T n+m-2 C(^ 1,171 1))T ••• • 


(4.7) 


The difference in comparison to the first moment is that, by the arguments above, the 
additional term ((n— 1, m— 1) is of the same order in N as the supposed leading term. 
This is also in contrast to the outcome for the second order term in the eigenvalue 
representation (14.31) . where the effect of removing the matrix diagonal leaves only a 
1 /N correction. Nevertheless we present arguments that allow for it to be neglected. 
Let us continue by inserting the expressions (14.51) . (14.G|) and (14.71) into the appropriately 
scaled version of the n simultaneous equations m , which determine the stationary 
solution Q (the method for calculating the error terms in the analogous eigenvalue 
representation approach is discussed at length in m and thus we refrain from details 
here). Therefore, for large N, the stationary solution Q for approximately satisfies 



To estimate the contribution of£(n — l,m — 1) we replace the exact value with its 
mean, i.e. TV -1 J2 p Bp^B™- 1 ss r n _ir m _i. For all matrices B £ IB n we have t\ = 0 
and 72 = N(N — 1)/N 2 = 1 — 1/N, meaning our space of variables is reduced to r n 
for n = 3,... N. Assuming then, that in all the remaining directions our JPDF Q is 
constant (as is the case in the GOE expression (12.51) 1 we find for n > 3 



where we have used that 9r m _i/9r m = 0 and dr n _i/dr m = <5 mi n_i for all n,m < N. 
This results in a term which is of order 1/N less than the corresponding term on the 
LHS and a full order 1 /N 2 less than the leading term. 

5. Discussion 

The efforts invested in developing the formalism presented above were motivated by 
our initial observations regarding random regular graphs. Dyson’s original model 
could not be transcribed to this matrix ensemble as the perturbation formula is 
effectively useless (a consequence of small separation between eigenvalues) in this 
context (see [24|1. Here we offer a method which does away with the requirement of 
the perturbation formula and therefore offers a potential method for circumventing 
such problems. We have demonstrated this method in the standard Gaussian setting 
and also illustrated how this can be used for Bernoulli matrices. The former case leads 
immediately to two previously unseen identities regarding symmetric functions, which 
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are proved directly below. Finally we also note the relation with those studies [M3] 
regarding the distributions of traces. Except for in. these works did not consider 
any dynamical aspects and so what we have outlined here may offer alternative ways 
for studying traces distributions. For instance one should be able to apply the same 
techniques to the circular ensembles. 
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Appendix A. Proof of Proposition [T] by Don Zagier 

Following the notation of the paper, we let X a (a = 1,..., N ) be independent variables 
and let Cj (0 < i < N ), t n (n = 0,1,...) and A (discriminant) be the elements of the 
algebra S = C[Ai,..., Avr] Siv of symmetric polynomials in the A a defined by 

N N 

*(x) ■■= n (*-*«) = E c * xz ’ 

ol— 1 i —0 

N 

t n = £ A" , A = disc($) = (A a — \p) 2 . 

ct=l 1 <cx.</3<N 

For n < 0 we set t n = 0. We have Cn = 1 and t 0 = N, while both (ci,..., cjv) and 
(ti,..., tjv) generate the algebra S. In particular, if we take the latter as coordinates 
on S, then we can ask for the values of dt n /dt m and dA/dt m for n > 0 and 1 < m < N. 
(Of course the former is 5 nrn for 0 < n < N, so it is only interesting if n > N.) 

The identities (EH) and m were proved in the body of this paper using an 
indirect proof coming from random matrix theory. Here we give a purely algebraic 
verification of both of these identities, and some small generalizations. For the reader’s 
convenience we repeat these identities here, expressing the second one in terms of the 
polynomial invariant A rather than its square-root G. 

Proposition. For n > 0 we have 

= y] titj + (n+l)t„, (A.l) 

i,i> 0 
i+j=n 

= ^2 titj ~ ( n + 1) t„ ■ (A.2) 

i,j> 0 

i+j=n 

We use that the logarithmic derivative of 4>(X) is a generating series for the £ n , 


N 


m—1 
N 


dt 


m ■ 


n+m 


8t r 


1 UA <9A 

A ^ mtn+m dtZ 

m=l 
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where the last expression can be taken either as a formal power series in S^l/X]] or as 
a holomorphic function in the annulus |X| > max Q |A a | if the A a are complex numbers. 
Dividing (IA.1I) and (IA.2I) by X" +2 and summing over n > — m (or equivalently n > 0, 
since dt n + m /dt m vanishes for —m < n < 0), we can rewrite these two identities as 

2 E m d T( x ) X m-1 = T ( X }2 _ (A.3) 

\ Otm 

m= 1 

and 

EE V = T(X) 2 + T'(X). (A.4) 

A ^ 9t m 

m= 1 

For the proof, we define polynomials <J>q,(X) and coefficients c a , n for 1 < a < N and 
0 < n < N — lby 




«(*) = ] 




X — 1 <F(X) 

A q -A^ _ $'(A a ) X - A q 


jv-i 

E< 

n —0 


\rr 


Then < h a (A / g) = S ai g, so {c a ^ n ) is the inverse of the Vandermonde matrix (A^)„ ia - On 
the other hand, we have A = A™ -1 , so c Q , m _i = m for 1 < m < N. Hence 


dT{X) ^ dT{X) _ ^ Ca,m—1 

dt m 1 dXa E (X — Aq,) 2 


(A.5) 


so each term dT(X)/dt m is a rational function of the form P m (X)/<f>(X) 2 
where P m (X) is a polynomial of degree < 2n — 2. Multiplying (IA.5I) by X m_1 and 
summing over m = 1,,N gives 


N 


E 


d T(X) 

dt m 


I" 1 " 1 


y- 3*a(X) 

E (X-A Q ) 2 

a—1 x 

N ^ 1 1 

(a-a q )3 


N 

*(X) J2 Res z _A„ 

a—1 


1 dz A 
$(I) (X - zf ) 


*(X) Re s ,- X (j—X_) 

$(X) d 2 1 1 $"(X) 

2 dX 2 $(X) “ _ 2 <F(X) 
T'(X) , T(X) 2 
2 + 2 ’ 


$'(X) 2 

4>(X) 2 


where in the fourth line we have used the residue theorem. This prove the first 
identity (IA.3I) . The calculation for (IA.4I) is similar. We have 


m dA 
2A di~ 


1 

2 


N 


a=l 


d log A 
d\ a 


N 


Xa _ x 

ppot 


a—1 


N 


^ ^ ^a,m— 1 
a=l 


$a(Aq) 

^a(Aa) ’ 
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Substituting into this the identity 


^(Aq) 

$a(A a ) 


m 


i 


t-Xn 


t=A, 


f 4>'(A Q + e) 1 

V ^(Aa + e) 


< h , (A Q ) + $"(A q ) e + ••• 1 

4>'(A a )e + A$"(A a )e 2 + ••• e 


£=0 


£ — 0 

1 d>"(A a ) 

2 $'(A Q ) 


multiplying by X m 1 and summing over m, we obtain the second identity (IA.4I) : 


1 y' m dA X m-i _ y' $ "(Aq) ^ cyI - V 

a \ m a^ x ~ 


N 


= $(X) Res z=> 

a=l 

= $(X) Res, =x 
$"(X) 


$"(z) dz 

<b'(z) (X - z)$(z) 

$"( 2 ) 


$'(X) 


= T(X) + 


(z - X)$(z)$’(z) 
T'{X) 


T{X ) ■ 


□ 


We mention that one can use the same method of calculation to obtain other 
identities of this type. For instance, 

JV N 




m=l 


a=l 

N 


K(X) 

(x - a q ) 2 


= E*7 


$'(X) 


4>(X) 


y^a(Aa) V(X-A a ) 3 (X — A a ) 4 
$'(X) $(X) \ dz 


= R.es z= * 


(i 


0 : 


V(z-x ) 3 (z-x) 4 y$( z )_ 

&(x) ( 1 V' , / 1 \ 

2 V$(X)/ 6 V$(X)/ 


,4>(X) 

= — t(x ) 3 - ir"(x) 


6 


and hence, in analogy with (I A. 1 1> . 

N 

3 'y^ m{m — 1) 


dt m 


+ + + (n + l)(n + 2) + 

— / J LiLjZk Z n 


i, j, k> 0: 
i+j + k = n 


Identities with polynomials of higher degree in m on the left could be proved in the 
same way. 
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